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Two numerical algorithms for the computation of eigenvalues of Dirac operators in lattice gauge theories are 
described: one is an accelerated conjugate gradient method, the other one a standard Lanczos method. Results 
obtained by Cullum's and Willoughby's variant of the Lanczos method (whose convergence behaviour is closely 
linked with the local spectral density) are presented for euclidean Wilson fermions in quenched and unquenched 
SU(2) gauge fields. Complete spectra are determined on lattices up to 8'^ • 12, and we derive numerical values for 
fermionic determinants and results for spectral densities. 



1. INTRODUCTION 

In order to study questions of chiral symme- 
try breaking [1-6] and universality [7,8], and also 
in the context of Liischer's fermion algorithm [9- 
11], one is interested in the eigenvalues of (75 
times) the gauge covariant Dirac operator which 
are close to the origin. In this talk we focus on the 
numerical aspects of the determination of spectral 
properties of lattice Dirac operators. Two ver- 
sions of lattice fermions are in use today [12]: Wil- 
son fermions [13] and staggered fermions [14,15]. 

Low-lying eigenvalues of staggered fermions 
were investigated e.g. in [2]. The present au- 
thor obtained complete spectra [16], and these 
data were further analysed by Halasz and Ver- 
baarschot [8]. They were interested in universal 
fluctuations in spectra of lattice Dirac operators, 
and they verified that the correlations among the 
eigenvalues of staggered fermions in SU(2) gauge 
fields are described by the Gaussian symplectic 
ensemble of random matrix theory with the chi- 
ral symmetry of the Dirac operator built in. 

Eigenvalues of Wilson fermions were computed 
e.g. in Refs. [17,3], and recently in [18]. Some 
of the latter results will be presented in this talk 
where we consider the hermitean operator 

= 75(i? + m)/(8 + m) (1) 
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with periodic boundary conditions. Q is nor- 
malized such that its eigenvalues are between -1 
and 1. (_D -|- m) is the Dirac operator for Wilson 
fermions of bare mass m. On a four-dimensional 
lattice A of sites x (with lattice spacing a = 1), 
[D-\-m) acts on a lattice spinor ^ as follows, "'^ see 
e.g. [12], 

Zk 

~2 ^^^^ ~ ^(*' * + ^^)<^b i^'^'^ix + n) 

+ (]1 + 7m)«/3 U(x,x- ii)ab i^^''(x - //)}. (2) 

Here k = (2m-|-8)~^ denotes the hopping param- 
eter and X ± /J, is the nearest neighbour site of x 
in ±//-direction. The gauge field U{x,x ± /j,) G 
SU(Nc) lives on the links (x, x ± /j,) of the lattice 
and is generated by some Monte Carlo process 
[12]. On the rhs of eq. (2) an implicit summation 
over the spinor indices {j3 = 1, . . . , 4) and colour 
(& = 1, . . . , A^c) is understood. 

Two algorithms for eigenvalues will be de- 
scribed in Sec. 2, and results in four-dimensional 
SU(2) gauge fields will be presented for complete 
spectra of Wilson fermions in Sec. 3. 



"^The hermitean euclidean 7-matrices 7^, /i = 1,2,3,4, 
satisfy the Clifford algebra {7^,71^} = 25^1^11. Moreover, 
75 = 71727374 anticommutes with all of them and 7^ = H 
is the 4x4 unit matrix. 
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2. ALGORITHMS FOR EIGENVALUES 

Maybe the most obvious candidate method 
which comes to mind to numerically determine 
eigenvalues of very large (hermitean or anti- 
hermitean) matrices^ is a Lanczos algorithm [19] 
which we explain in Sec. 2.2. Variants of this 
method have been used in lattice field theory for 
a long time, see e.g. Refs. [20,2] for staggered fer- 
mions and Refs. [17,3] for Wilson fermions. 

In order to obtain an overview of complete 
spectra the Lanczos method is well-suited in lat- 
tice QCD, as we will see in Sec. 3. However, 
the method can be problematic whenever one is 
interested in only a few eigenvalues. Whether 
eigenvalues have converged can be estimated only 
from experience [18], and not from a rigorous er- 
ror bound. ^ Moreover, a Lanczos method cannot 
provide information about multiplicities. 

Some applications also require knowledge of the 
eigenvectors, for instance to isolate the contri- 
bution of low-lying eigenmodes to physical ob- 
servables. In such cases a conjugate gradient 
(CG) algorithm is superior, and even if one just 
wants knowledge of eigenvalues the CG method 
is favourable because it does not suffer from the 
drawbacks and potential pitfalls of the Lanczos 
method. So let us turn first to an accelerated 
version of a straightforward CG implementation 
[22]. 

2.1. Accelerated CG method 

In [23] a CG method was proposed for the com- 
putation of low-lying eigenvalues of Q^, cf. also 
[24,25]. More generally, this procedure can be 
used to obtain extreme eigenvalues of a hermitean 
operator. The CG method is based on the gauge 
covariant extremization of the Ritz functional'* 



(^^0) 



(3) 



where for the k-th lowest eigenvalue ip is kept 
orthogonal to the eigenspace of the (k — 1) 
lower eigenvalues. The pure CG minimization of 
Ref. [23] was accelerated in [22] by alternating 



is represented by an n X n matrix with n = 4A''c|A|. 
■^In principle one can quote error bounds [19,21], but they 
are not practical in lattice gauge theories. 
^ ( • , • ) denotes the scalar product of the Hilbert space. 



CG searches with exact diagonalizations in the 
subspace spanned by the numerically computed 
eigenvectors. This method is attractive because 
of the following key features. 

• Rigorous error bounds for eigenvalues can 
be derived just from the last CG iterate. 

• The correct multiplicities are detected. 

• Approximations to eigenvectors are ob- 
tained as a by-product. 

• The pure CG algorithm is speeded up 
through the intermediate diagonalizations 
by a factor of 4 — 8. 

Furthermore, the algorithm is numerically very 
stable, even if one is restricted to single precision 
arithmetics, and it can be parallelized in an effi- 
cient way. For more details we refer to Ref. [22]. 

2.2. Lanczos method 

The Lanczos procedure is a technique that can 
be used to solve large, sparse, symmetric or her- 
mitean eigenproblems® [19]. The idea is to trans- 
form a given hermitean nxn matrix A into a sim- 
ilar symmetric tridiagonal matrix T = V~^AV 
with unitary V , and then T is diagonalized. The 
transformation of A can be performed iteratively. 
If one writes V = [vi,V2, ■ ■ ■ ,Vn) with column 
vectors Vi ("Lanczos vectors") and 



T 



( a\ /?! 

V 



/3n-l 

an 



(4) 



/3n-l 

then AV = VT is equivalent to (i = 2, 

Avi = aivi -I- (3iV'2 , 



Avi 
Av„ 



(5) 



Given an initial — generally random — vector vi 
and using the orthonormality among the Vi with 
respect to the canonical scalar product ( • , • ), one 
can determine iteratively from these equations: 



^In Refs. [20,17,3] a non-hermitean Lanczos method was 
used and it was concluded that it works on small lattices. 
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Ui = {vi,Avi), I3'f = {{vi,A'^Vi)-al-(3f_^), (with 
/3o = 0), and Vi+i = (3~^{Avi - Pi^iVi^i - UiVi) 
as long as j3i ^ (otherwise the iteration is 
stopped). T will be a real matrix also in case 
that A is complex, and the sign of j3i is arbitrary. 

In exact arithmetic, the Lanczos iteration 
should finish after at most n steps and the last 
equation in (5) would be automatically fulfilled. 
For this case there exists a convergence the- 
ory [19]. In practice, however, there are severe 
problems with a straightforward implementation 
[19,21]. These problems are caused by rounding 
errors and loss of orthogonality among the Lanc- 
zos vectors.® In principle the latter problem can 
be circumvented by storing all the Lanczos vec- 
tors and enforcing orthogonality among them by 
hand. But then one is restricted to small lattices 
because of computer memory or I/O limitations. 

A practical Lanczos method is due to Cullum 
and Willoughby [21]. In their proposal one per- 
forms no reorthogonalization, and one continues 
the iteration (5) for an a priori unspecified count. 
In this way a sequence of j x j tridiagonal matri- 
ces T'--'-', j = 1,2,... is generated. There exists 
a recipe [21] how the eigenvalues of T'--'-' are con- 
nected with eigenvalues of Q, see [18] for details. 

3. NUMERICAL RESULTS 

Cullum's and Willoughby's Lanczos method 
[21] was studied in quenched and unquenched^ 
SU(2) gauge fields, i.e. #c = 2 in the following. 
Tridiagonal matrices were diagonalized by means 
of the Numerical Recipes routine "tqli" [26] which 
implements the QL algorithm with implicit shifts. 
Two eigenvalues A were counted as different when 
they differed by more than 10"^*^. This number 
is arbitrary but it is chosen such that it is small 
compared with the gaps in the spectra, and large 
compared with round-off errors. The computer 
program was checked for gauge covariance, and 
it was also verified that free spectra are obtained 
correctly, except for multiplicities. 

In the following we present results for complete 
spectra. Further results and a comparison of the 
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Figure 1. Convergence of the Lanczos method. 

partially converged Lanczos method with the ac- 
celerated CG algorithm for low-lying eigenvalues 
of Ref. [22] can be found in [18]. 

3.1. Convergence behaviour 

Let us start by looking at the convergence be- 
haviour of the Lanczos algorithm. We would like 
to concretize what is meant by the general state- 
ment of Ref. [21] that the local gap structure plays 
a role. An example with n = 20 736 eigenvalues 
in an unquenched configuration on a 6^ • 12 lat- 
tice at® l3 = 2.12, K = 0.15 is shown in Fig. 1. 
The Lanczos method was run both with A = Q 
and with A = . Every 500 iterations conver- 
gence was checked by a comparison with refer- 
ence data from T^S2944)^ p^^g 

ffive the number 



^This loss of orthogonality is not necessarily due to the 
accumulation of round-off errors [19]. 
^With two flavours of dynamical fermions. 



^Here and in the following /5 denotes the coupling constant 
of the gauge part of the action, see e.g. [12]. 
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j of Lanczos iterations in multiples of n which 
are required for convergence of A or A^, and di- 
amonds indicate 1.75[Aniax — -^min — l'^|]p('^) ^nd 
3.5[(Amax-Amin)A-A2]p(A2), where Xmax/mm de- 
notes the highest /lowest eigenvalue of Q. 

We notice that the number of Lanczos itera- 
tions required for convergence is intimately re- 
lated with the local spectral density. For this 
reason it is advantageous to diagonalize the un- 
squared operator in case that one requires the 
lowest or highest eigenvalues of . This state- 
ment will become more and more significant as a 
critical point is approached in the /3-K-plane. 

3.2. Complete spectra 

Complete spectra were determined in the con- 
figurations quoted in Tables 1 and 2 below. 
In nontrivial gauge fields no degeneracies were 
found, neither for Q nor for . In case of Q 
we always encountered an equal number of neg- 
ative and positive eigenvalues. Furthermore, no 
discrepancies were found when the squared eigen- 
values of Q were compared with the results for . 

Examples for complete spectra in a quenched 
gauge field at /3 = 1.8 on a G'* lattice are given 
in Fig. 2. The integrated densities of eigenvalues 
N{X) and N{X^) follow directly from the numer- 
ical data. They are normalized such that they 
take values between zero and one. At k = 0.15 
curves of a quenched 12'* gauge field^ at /3 = 1.8 
and unquenched gauge fields at /3 = 2.12 on 6^-12 
and 8^ • 12 lattices coincide (on the scale of the 
figure) with the result on the G'* lattice. 

We have the following consistency checks which 
provide good evidence that all computed com- 
plete spectra are correct. First, on all investi- 
gated 4^* - 8^ • 12 lattices the correct number of n 
different eigenvalues was found. Second, we have 
analytical sum rules for powers of the eigenval- 
ues, cf. also Ref. [16]. The traces (over colour 
and spinor indices) of Q and equal zero, and 
the trace of reads in any unitary gauge field 

TrQ2 = 4N,|A|(4+-L)/(4+^)2 . (6) 

We obtained |TrQ|<10-8, |TrQ3|<10-^, and 

Because of memory limitations, on the 12* lattice all but 
four of the 165888 eigenvalues were found; see [18]. 
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Figure 2. Integrated spectral densities. 

Tr came out with a relative accuracy of 10~® - 
10~^^ (decreasing with increasing |A|). One could 
check for further sum rules by examining traces 
of higher powers of Q. Finally, we confirmed that 
the effort for the determination of complete spec- 
tra grows with |Ap, j = An worked in all cases. 

3.3. Fermionic determinants 

For free Wilson fermions spectra and determi- 
nants are known analytically [18]. Quenched / 
unquenched examples are contained in Table 1/2. 
One concludes a nice exponential dependence of 
the determinant on the lattice volume, already for 
relatively small lattices. 

In a quenched Monte Carlo simulation the 
fermionic determinant is kept at a fixed value. 
From Table 2 one can get a feeling for the fiuc- 
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Table 1 

Examples for Wilson fermions in quenched SU(2) gauge fields: (logj^g det Q) / {4:Nc\A\) . 



|A| 


/9 


K = 0.1 


K = 0.125 


K = 0.15 


K = 1/6 


K = 0.20 


K = 0.25 


6* 


1.80 


-0.254757 


-0.299665 


-0.339289 


-0.362812 


-0.402336 


-0.440477 


6* 


2.80 


-0.254331 


-0.298278 


-0.334999 


-0.354646 


-0.383828 


-0.409753 


6* 


0.00 


-0.255279 


-0.301046 


-0.342455 


-0.368026 


-0.415065 


-0.477151 


8* 


0.00 


-0.255277 


-0.301041 


-0.342445 


-0.368010 


-0.415038 


-0.476955 



Table 2 

Examples for Wilson fermions in unquenched 
SU(2) gauge fields: (log^o det 0)/(4#c|A|). 



p{X), /3=1.8 



tuations of the determinant in an unquenched 
run. We can compare the values of eight indepen- 
dent dynamical 6^-12 configurations at /3 = 2.12, 
K = 0.15. The logarithmic entries in Table 2 fiuc- 
tuate by « 0.0003 which translates to a fiuctu- 
ation of the determinant itself by six orders of 
magnitude on the 6^ • 12 lattice. 

3.4. Spectral densities 

The density of eigenvalues of around zero 
can be related with the chiral limit, if one con- 
nects the spectral density of with that of the 
Dirac operator (not multiplied by 75) in the spirit 
of the "pion norm" of Ref. [1]. Such an analysis 
will however be done elsewhere. From the algo- 
rithmic point of view spectral densities p{X) of 
Q and p{X^) of are interesting because they 
determine convergence properties (cf. Fig. 1). 

Figs. 3 and 4 show normalized quenched p's 
on a G'* lattice at /3 = 1.8 and in a random S'* 
gauge field. (Additional results at /3 = 2.8 are in 
[18].) The solid curves indicate reference data of 
unquenched gauge fields at k = 0.15, /3 = 2.12 on 
6^ • 12 and 8^ • 12 lattices, where results coincide. 



1.2 



|A| 


/9 


K = 0.15 




4* 


1.75 


-0.338947 


0.8 


6^ • 12 


2.12 


-0.337169 


6^ • 12 


2.12 


-0.337129 


0.6 


6^ • 12 


2.12 


-0.337282 




6^ • 12 


2.12 


-0.337388 


0.4 


6^ • 12 


2.12 


-0.337385 




6^ • 12 


2.12 


-0.337266 


0.2 


6^ • 12 


2.12 


-0.337269 




6^ • 12 


2.12 


-0.337235 





8^ • 12 


2.12 


-0.337295 






p(A2), /3=1.8 
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Figure 3. Spectral densities at /3 = 1.8. 

The information in Figs. 3 and 4 is redundant 
because of the relation p{X^) = p{X)/{2X), which 
can be verified numerically to a very high pre- 
cision. Despite this redundancy we consider it 
worthwhile to present results for p{X) as well as 
for p{X^). Note for instance that p{X^) stays fi- 
nite (must diverge) as — when p{X = 0) = 
{p{X = 0) is finite). 
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p{X), /3 = 0.0 




Figure 4. Spectral densities at (] 



I wish to thank B. Bunk, K. Jansen, H. 
Simma, and U. Wolff for discussions. C. Liu 
and K. Jansen are thanked for providing the un- 
quenched Wilson configurations. Financial sup- 
port by Deutsche Forschungsgemeinschaft under 
grant Wo 389/3-1 is gratefully acknowledged. 
The computations reported here were performed 
on HP workstations of the Humboldt-University. 

REFERENCES 

1. K.M. Bitar, A.D. Kennedy, and P. Rossi, 
Phys. Lett. B234 (1990) 333. 

2. S.J. Hands and M. Teper, Nucl. Phys. B347 
(1990) 819. 

3. L Barbour, E. Laermann, T. Lippert, and K. 
Schilling, Phys. Rev. D46 (1992) 3618. 



7. 
8. 

9. 

10 



11 



12 



13 



14 

15 

16 

17 

18 
19 



20 



21 

22 

23 

24 
25 
26 



H. Leutwyler and A. Smilga, Phys. Rev. D46 
(1992) 5607; A.V. Smilga, hep-th/9503049. 

A. Hoferichter, V.K. Mitrjushkin, and M. 
Miiller-Preussker, hep-lat/9506006. 

F. X. Lee, H.D.Trottier, and R.M. Woloshyn, 
hep-lat/9509028. 

J.J.M. Verbaarschot, hep-th/9405006. 

M.A. Halasz and J.J.M. Verbaarschot, Phys. 

Rev. Lett. 74 (1995) 3920. 

M. Liischer, Nucl. Phys. B418 (1994) 637. 

B. Bunk, K. Jansen, B. Jegerlehner, M. 
Liischer, H. Simma, and R. Sommer, Nucl. 
Phys. B (Proc. Suppl.) 42 (1995) 49. 

C. Alexandrou, A. Borelli, Ph. de Forcrand, 

A. Galli, and F. Jegerlehner, hep-lat/9506001. 
L Montvay and G. Miinster, Quantum fields 
on a lattice (CUP, Cambridge, 1994). 

K.G. Wilson, Erice 1975, in: New phenomena 
in subnuclear physics, ed. A. Zichichi (Plenum 
Press, New York, 1977). 

J. Kogut and L. Susskind, Phys. Rev. Dll 
(1975) 395. 

H. Sharatchandra, H. Thun, and P. Weisz, 
Nucl. Phys. B192 (1981) 205. 
T. Kalkreuter, Phys. Rev. D51 (1995) 1305; 
Nucl. Phys. B (Proc. Suppl.) 42 (1995) 882. 
R. Setoodeh, C.T.H. Davies, and LM. Bar- 
bour, Phys. Lett. B213 (1988) 195. 
T. Kalkreuter, hep-lat/9509071. 

G. H. Golub and C.F. v. Loan, Matrix compu- 
tations, second edition, (The Johns Hopkins 
University Press, Baltimore, 1990). 

LM. Barbour, N.-E. Behilil, P.E. Gibbs, G. 
Schierholz, and M. Teper, in: The recursion 
method and its applications. Springer Series 
in Solid-State Sciences 58, eds. D. G. Pettifor 
and D. L. Weaire (Springer, Berlin, 1985). 
J. Cullum and R.A. Willoughby, J. Comp. 
Phys. 44 (1981) 329. 

T. Kalkreuter and H. Simma, hep-lat/- 
9507023, to appear in Comp. Phys. Comm. 

B. Bunk, K. Jansen, M. Liischer, and H. 
Simma, internal DESY-report (Sep. 1994). 
M. Geradin, J. Sound Vib. 19 (1971) 319. 

L Fried, J. Sound Vib. 20 (1972) 333. 

W.H. Press et al.. Numerical recipes (CUP, 

Cambridge, 1992). 



